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1  Abstract 

The  objective  of  this  project  is  to  develop  and  analyze  stable  time  discretizations  suitable  for  the 
simulation  of  hyperbolic  time-dependent  partial  differential  equations.  Implicit  and  explicit  multi- 
step  multi-stage  time  discretizations  with  optimal  time-step  restrictions  have  been  developed,  as  well 
as  implicit  Runge-Kutta  methods  with  downwinding  and  unconditional  stability.  New  directions 
have  been  explored  for  implicit-explicit  methods,  which  are  useful  for  problems  with  convection 
and  diffusion.  Finally,  novel  provably  stable  multi-step  time  discretizations  for  use  with  Fourier 
pseudo-spectral  spatial  approximations  of  the  three  dimensional  viscous  Burgers  equations  and 
Navier-Stokes  equations  have  been  developed. 

2  Summary 

In  this  second  year  of  the  grant,  I  have  made  progress  in  three  directions.  First,  we  developed 
optimization  codes  that  search  for  implicit  and  explicit  multi-stage  multistep  methods  of  more  than 
two  steps.  These  codes  are  now  being  utilized  to  find  three  step  Runge-Kutta  methods  with  optimal 
SSP  coefficients.  Next,  we  have  made  significant  progress  in  developing  Runge-Kutta  methods  with 
downwinding  and  unlimited  time-step,  as  well  as  testing  appropriate  boundary  conditions  for  such 
methods.  We  investigated  implicit-explicit  methods  with  interesting  SSP  properties,  a  promising 
direction  we  intend  to  continue  to  develop.  Finally,  we  have  developed  a  series  of  multistep  methods 
which  are  provably  stable  for  use  with  Fourier  pseudospectral  spatial  discretizations  of  3D  viscous 
Burgers’  equation  and  the  3D  Navier-Stokes  equations.  Also  in  this  year,  we  published  a  book  on 
SSP  time  discretization  methods. 

3  Objectives  and  accomplishments 

1.  Objective:  Study  the  class  of  high  order  implicit  and  explicit  SSP  Runge-Kutta 
methods  with  downwinding  to  allow  for  higher  order  SSP  methods. 

Motivation:  To  more  easily  analyze  SSP  methods,  we  rewrite  Runge-Kutta  methods  in  the 
form: 

U0)  =  un, 

l—l 

u(i)  =  ^(^uW  +  Atft^uW)),  ai)k  >  0,  i  =  (3.1) 

k= 0 
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un+ 1  =  u(m). 

Explicit  SSP  Runge-Kutta  methods  are  known  to  be  limited  to  fourth  order  and  implicit  SSP 
Runge-Kutta  methods  are  limited  to  sixth  order.  However,  if  we  allow  the  use  of  negative  co¬ 
efficients  Pi  k  it  is  possible  to  overcome  this  order  barrier.  The  presence  of  negative  coefficients 
requires  the  use  of  a  modified  spatial  discretization  for  these  instances.  When  fii.k  is  negative, 
Pi,kF(u^)  is  replaced  by  PitkF(u^),  where  F  approximates  the  same  spatial  derivative(s)  as 
F.  but  the  strong  stability  property  holds  for  the  first  order  Euler  scheme,  solved  backward 
in  time.  Numerically,  the  only  difference  is  the  change  of  the  upwind  direction. 

The  order  barrier  on  explicit  SSP  Runge-Kutta  methods,  and  the  bounds  on  the  SSP  coeffi¬ 
cient  for  implicit  SSP  Runge-Kutta  methods,  leave  no  efficient  options  for  methods  of  order 
five  and  above.  This  gap  may  be  filled  by  explicit  and  implicit  SSP  Runge-Kutta  methods 
with  downwinding. 

Accomplishments: 

(a)  I  have  created  an  optimization  code  in  MATLAB  which  seeks  implicit  methods  with 
downwinding  with  a  large  allowable  SSP  coefficient. 

(b)  The  code  has  been  re-written  in  Python  to  facilitate  open  source  reproduction  of  results 

(c)  This  code  has  produced  results  for  implicit  methods  of  orders  2-4,  and  in  future  will 
produce  higher  order  methods 

(d)  A  similar  optimization  code  has  been  developed  for  explicit  methods  and  for  diagonally 
implicit  methods. 

(e)  These  methods  have  been  tested  on  a  test  suite  of  specially  selected  of  problems. 

2.  Objective:  Study  classes  of  high  order  implicit  and  explicit  general  linear  methods 
to  find  SSP  methods  with  optimal  time-step  restrictions. 

Motivation:  Without  the  use  of  downwinding ,  explicit  SSP  Runge-Kutta  methods  are  limited 
to  fourth  order  and  implicit  SSP  Runge-Kutta  methods  are  limited  to  sixth  order.  SSP  multi- 
step  methods  do  not  suffer  from  this  order  barrier,  but  have  very  restrictive  SSP  coefficients. 
Efficient  explicit  SSP  methods  of  order  greater  than  four  are  frequently  desirable,  particularly 
when  dealing  with  high  order  spatial  discretizations.  General  linear  methods,  which  have 
multiple  steps  and  multiple  stages  have  the  potential  to  combine  the  properties  of  multistep 
and  Runge-Kutta  methods,  and  so  provide  an  advantage  over  these  methods  by  allowing  a 
larger  step-size  [4,  3].  We  have  shown  [1]  that  explicit  general  linear  methods  have  a  bound 
on  the  SSP  coefficient  which  is  equal  to  the  number  of  stages.  Even  considering  this  bound, 
explicit  general  linear  methods  may  be  found  that  have  order  p  >  4  and  larger  SSP  coefficient 
than  the  multistep  methods.  A  first  effort  in  this  direction  has  been  made  by  [2],  which 
considered  a  subset  of  general  linear  methods. 

Multistep  Runge-Kutta  methods  are  a  straightforward  generalization  of  Runge-Kutta  and 
linear  multistep  methods,  and  take  the  form 


yni  =  Ej=i  dijUn+1-j  +  AtJ2j=i  !<*<«, 

un+l  =  Ef= 1  j  +  A  t  v  •  ,  hjfuj'j). 

Here  the  values  un  denote  solution  values  at  the  times  t  =  nAt,  while  the  values  y"  are 
intermediate  stages  used  to  compute  the  next  solution  value.  We  will  also  consider  a  simple 
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generalization  of  these  methods,  based  on  the  following  reasoning.  For  some  methods,  it  may 
happen  that  the  row  i  of  A  is  identically  zero  and  row  i  of  D  is  (1,  0, . . . ,  0),  so  that  y™  =  u11. 
Then  the  method  involves  f{un ),  and  at  any  step  we  will  have  computed  already  f(un+1~^) 
for  j  =  1, . . . ,  k,  so  these  values  may  as  well  be  used  in  computing  the  next  step.  This  leads 
to  methods  of  the  form 

Vi  =  un, 

yni  =  EjLi  dijU"+l-j  +  A t  E -=2  Oijfivr+'-i)  +  At  E -=i  o^/fo?),  2  <  i  <  s, 

Un+1  =  EjU  ejUn+1~i  +  At  E*=2  bjf(un+1~j )  +  At  E*=1  bjf(y?). 

This  form  is  more  suitable  for  finding  explicit  methods. 

Accomplishments: 

(a)  We  have  found  high  order  (up  to  9th)  explicit  SSP  RK  methods  with  large  effective  SSP 
coefficients  among  the  class  of  2-step  RK  methods.  The  paper  has  been  published. 

(b)  We  currently  have  a  code  running  to  numerically  find  optimal  high  order  implicit  SSP 
multistep  RK  methods  with  more  than  2  steps.  We  have  verified  numerically  that  the 
optimal  methods  are  all  one-step  diagonally  implicit  methods  with  c  =  2.  The  paper  is 
in  preparation. 

(c)  We  modified  the  above  code  running  for  diagonally  implicit  SSP  multistep  RK  methods 
with  more  than  2  steps,  and  will  be  running  it  next. 

(d)  We  modified  the  above  code  running  for  explicit  SSP  multistep  RK  methods  with  more 
than  2  steps  using  second  form  above.  The  code  has  produced  optimal  methods  up  to 
order  6  and  is  currently  working  on  methods  of  orders  7-10.  The  paper  is  in  preparation. 

3.  Objective:  Study  possibilities  for  efficient  implicit-explicit  SSP  methods.  Fre¬ 
quently  we  are  faced  with  problems  of  the  type 

un+l  =un  +  AfF(un )  +  A*G(U"+1) 

which  are  unconditionally  stable,  but  only  first  order.  Because  higher  order  implicit  methods 
(without  downwinding)  cannot  be  unconditionally  unstable,  we  cannot  extend  this  to  higher 
order. 

However,  if  we  are  given  a  different  type  of  method,  one  for  which 

un+l  =  un  +  +  A^G^n+ 1) 

is  provably  unconditionally  stable,  then  we  can  develop  SSP  methods  which  are  higher  order 
and  unconditionally  SSP. 

Accomplishments: 

(a)  The  observation  above  is  an  accomplishment  -  it  changes  the  usual  paradigm  for  looking 
at  IMEX  methods. 

(b)  We  have  applied  this  idea  to  problems  with  artificial  viscosity  in  spectral  methods 

4.  Objective:  Tailor  semi-implicit  methods  for  the  stable  solution  of  the  3D  viscous 
Burgers’  equation  and  the  Navier-Stokes  equations 
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Motivation:  We  consider  the  three-dimensional  viscous  Burgers’  equation  with  periodic  bound¬ 
ary  conditions, 

dtu  +  u- Vm  =  vAu,  (3-2) 

where  u  =  ( u ,  v,  w)T  is  the  velocity  field  in  the  x ,  y,  and  z  directions,  and  v  >  0  is  the  viscosity. 
We  discretize  the  spatial  terms  using  a  Fourier  collocation  (pseudospectral)  method,  and  the 
time-derivative  using  a  variety  of  suitably  chosen  semi-implicit  multistep  methods.  Although 
this  method  has  been  widely  applied  to  this  problem  (and  the  similar  problem  of  the  Navier- 
Stokes  equations),  a  fully-discrete  stability  analysis  has  not  been  performed.  Moreover,  it  is 
not  clear  what  time-stepping  is  stable. 

We  begin  with  a  first-order  in  time  method  which  treats  the  nonlinear  convection  term  ex¬ 
plicitly  for  the  sake  of  numerical  convenience,  and  the  diffusion  term  implicitly  to  avoid  a 
severe  time-step  restriction: 

n.n+l  _  1ln 

- = -  +  un  •  V  Nun  =  |/Ajv«n+1,  (3.3) 

where,  for  example,  the  first  component  of  the  nonlinear  convection  is 

un-VNun  =  unVNxun  +  vnVNyun  +  wnVNzun. 

In  our  recent  paper,  we  performed  a  fully-discrete  analysis  and  found  that  this  method  is 
unconditionally  stable  for  any  final  time  T*  provided  that  the  time  step  At  and  the  space 
grid  size  h  are  bounded  by  given  constants 

At  <  Li(T*,  v),  h  <  L-2 (T* ,  v) . 

These  convergence  constants  depend  on  the  exact  solution,  as  well  as  T*  and  v. 

If  desired,  we  could  use  a  fully  explicit  method  and  prove  stability  for  a  time-step  restriction 
of  the  form  At  <  C Ax2,  at  which  point  we  could  apply  an  SSP  method  up  to  fourth  order 
and  maintain  the  stability  of  this  method.  However,  if  we  wish  to  avoid  this  approach  we 
need  to  tailor  the  time-discretizations.  (This  is  due  to  the  problem  that  IMEX  methods  do 
not  maintain  the  SSP  property,  as  mentioned  above). 

Accomplishments: 

(a)  We  derived  a  provably  stable  second  order  scheme  for  the  3D  Burgers’  equation.  For 
the  convection  term  we  use  a  standard  second  order  Adams-Bashforth  extrapolation 
formula,  which  involves  the  numerical  solutions  at  time  node  points  fn,  tn_1,  with  well- 
known  coefficients  3/2  and  —1/2,  respectively.  The  diffusion  term  is  treated  implicitly, 
using  a  second  order  Adams-Moulton  interpolation.  Treating  the  diffusion  term  with  a 
standard  second  order  Adams-Moulton  formula  leads  to  a  Crank-Nicolson  scheme,  which 
gives  us  difficulties  in  the  stability  analysis. 

Instead,  we  look  for  an  Adams-Moulton  interpolation  such  that  the  diffusion  term  is  more 
focused  on  the  time  step  tn+1,  i.e.,  the  coefficient  at  time  step  tn+1  dominates  the  sum 
of  all  other  diffusion  coefficients.  We  discovered  that  the  Adams-Moulton  interpolation 
which  involves  the  time  node  points  fn+1  and  gives  the  corresponding  coefficients  as 
3/4,  1/4,  respectively,  which  satisfies  the  unconditional  stability  condition.  Therefore, 
we  formulate  the  fully  discrete  scheme: 

un+1~t  un  +  \un-VNun  -  l-un-l-VNun-1  =  uAn  (^un+1  +  \un -1)  .  (3.4) 

This  scheme  is  provably  stable  in  the  same  way  as  the  first  order  scheme. 
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(b)  We  derived  provably  stable  third  and  fourth  order  schemes  for  the  3D  Burgers’  equation. 
Using  similar  ideas  we  derived  third  and  fourth  order  in  time  schemes  for  (3.2).  The 
nonlinear  convection  term  is  updated  by  an  explicit  Adams-Bashforth  extrapolation 
formula,  with  the  time  node  points  tn,  in_1,  ...,  tn~k+l  involved  and  an  order  of  accuracy 
k.  The  diffusion  term  is  computed  by  an  implicit  Adams-Moulton  interpolation  with  the 
given  accuracy  order  in  time.  To  ensure  unconditional  numerical  stability  for  a  fixed 
time,  we  have  to  derive  an  Adams-Moulton  formula  such  that  the  coefficient  at  time 
step  tn+1  dominates  the  sum  of  the  other  diffusion  coefficients.  In  more  detail,  a  £:-th 
order  (in  time)  scheme  takes  the  form  of 


un+1  —  un 
A t 


fc-i 

+  YJBlun-VNun~i 

i= 0 


uA 


N 


(dqUu+1 


k- 1 

+  Y,D 

i= 0 


m 


u 


n-j(i ) 


(3.5) 


in  which  Bi  I^Tg1  are  the  standard  Adams-Bashforth  coefficients  with  extrapolation  points 
tn,  tn~ 1 ,  ...,  tn~k+1,  j(i)  \k~Q  are  a  set  of  (distinct)  indices  with  j(i)  >  0,  and  Dq, 
Dju\  l^1  correspond  to  the  Adams-Moulton  coefficients  to  achieve  the  k- th  order  ac¬ 
curacy.  Moreover,  a  necessary  condition  for  unconditional  numerical  stability  is  given 
by 


k- 1 


A)>£  I  Dj{i) 

i= 0 


(3.6) 


To  derive  an  Adams-Moulton  formula  for  the  diffusion  term,  whose  coefficients  satisfy 
the  condition  (3.6),  we  require  a  stretched  stencil.  In  particular,  for  the  third  order 
scheme,  it  can  be  shown  that  a  stencil  comprised  of  the  node  points  tn+ 1 ,  tn~l  and  fn~3 
is  adequate.  The  fully  discrete  scheme  can  be  formulated  as: 


un+l  —  ur 


At 


23  n  _  n  4 

-| - un-VNun - u 

12  3 


■VNun~l  +  -ur 


■V  Nur‘ 


=  uAn  (  -un+l  +  —un~l  -  — un~ 3 
1  3  12  12 


(3.7) 


For  the  fourth  order  scheme,  we  use  an  Adams-Moulton  interpolation  at  node  points 
tn+1,  tn_1,  and  fn-'  for  the  diffusion  term.  Combined  with  the  Adams-Bashforth 
extrapolation  for  the  nonlinear  convection  term,  the  scheme  is  given  by 


un+l  _  un 


At 


55  ^  „  59 

+  ~u  -VNu  --u 


Q7 

■Vnu71-1  +  -vf 


■VNun~2  -  -ur 


=  uA]y 


V 1152 


470  118 

■un  1  - 


1152 


1152 


un~5  + 


43 

1152 


-u 


n—7 


■VNur‘ 
(3.8) 


In  the  third  order  case,  we  observe  that  a  direct  application  of  the  Adams-Moulton 
formula  at  the  nodes  tn+1,  tn  and  tn~l  (corresponding  to  j (i)  =  i  in  the  general  form 
(3.5))  does  not  give  a  formula  with  the  stated  stability  property.  For  example,  a  “naive” 
combination  of  3rd  order  Adams-Bashforth  for  the  nonlinear  convection  and  Adams- 
Moulton  for  the  diffusion  term  results  in  the  following  scheme 


un+l  —  ur 


At 


23  n  n  4 

H - U  -V jyU - u 

12  3 


■vJVi*n_1  +  Y2ur 


■'VNur‘ 


5 

12 


=  vAn  (  -un+1  +  -un  -  ^ un~l 


1 

12 
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which  violates  the  stability  condition  (3.6).  Moreover,  numerical  experiments  also 
showed  that  this  method  is  in  fact  unstable.  This  case  highlights  the  need  to 
choose  an  appropriate  time-discretization  to  couple  with  the  pseudospectral  method. 


4  Dissemination 

We  continue  to  update  our  SSP  RK  web-site  to  disseminate  the  results  of  the  study.  This  site 
serves  as  an  online  catalog  of  all  the  methods  studied,  noting  which  are  most  successful,  and 
commenting  on  the  theoretical  properties  of  each,  and  on  which  performed  best  with  which  spatial 
approximation. 

In  Summer  2010,  I  organized  two  minisymposium  sessions  on  SSP  methods  at  the  SIAM  annual 
meeting  in  Pittsburgh. 

A  book  on  SSP  methods,  written  with  colleagues  David  Ketcheson  and  Chi- Wang  Shu  was 
published  by  World  Scientific  Press. 

Acknowledgment /Disclaimer  This  work  was  sponsored  (in  part)  by  the  Air  Force  Of¬ 
fice  of  Scientific  Research,  USAF,  under  grant/contract  number  FA9550-09-1-0208.  T  he  views  and 
conclusions  contained  herein  are  those  of  the  authors  and  should  not  be  interpreted  as  necessarily 
representing  the  official  policies  or  endorsements,  either  expressed  or  implied,  of  the  Air  Force  Office 
of  Scientific  Research  or  the  U.S.  Government. 


References 

[1]  S.  Gottlieb,  D.I.  Ketcheson  and  C.-W.  Shu.  High  Order  Strong  Stability  Preserving  Time 
Discretizations.  Journal  of  Scientific  Computing  ,  vol  38,  No.  3  (2009),  pp.  251-289. 

[2]  C.  Huang.  Strong  Stability  Preserving  Hybrid  Methods.  Applied  Numerical  Mathematics ,  doi: 
10.1016/j.apnum.2008.03.030,  2008. 

[3]  D.I.  Ketcheson,  C.B.  Macdonald  and  S.  Gottlieb.  Optimal  implicit  strong  stability  preserving 
Runge-Kutta  methods.  Applied  Numerical  Mathematics,  doi:  10.1016/j.apnum.2008.03.034. 

[4]  H.W.J.  Lenferink.  Contractivity-preserving  implicit  linear  multistep  methods.  Mathematics 
of  Computation,  56:177-199,  1991. 

Personnel  Supported  During  Duration  of  Grant 

Sigal  Gottlieb,  Professor  of  Mathematics,  UMass  Dartmouth. 

Saeja  Kim,  Associate  Professor  of  Mathematics,  UMass  Dartmouth. 

Zachary  Grant,  Undergraduate  Student,  UMass  Dartmouth. 

Sidafa  Conde,  Undergraduate  Student,  UMass  Dartmouth. 

Publications 

1.  S.  Gottlieb,  F.  Tone,  C.  Wang,  X.  Wang,  and  D.  Wirosoetisno,  Long  time  stability  of  a  classi¬ 
cal  efficient  scheme  for  two  dimensional  Navier-Stokes  equations.  SIAM  Journal  on  Numerical 
Analysis  (2012)  50,  pp.  126-150 

2.  D.I.  Ketcheson,  S.  Gottlieb,  and  C.  B.  Macdonald,  Strong  stability  preserving  two-step  Runge- 
Kutta  methods.  SIAM  Journal  on  Numerical  Analysis.  (2012)  49,  pp.  2618-2639. 


6 


3.  Sigal  Gottlieb  and  Cheng  Wang,  Stability  and  convergence  analysis  of  fully  discrete  Fourier 
collocation  spectral  method  for  3-D  viscous  Burgers’  Equation.  Submitted  to  Journal  of  Sci¬ 
entific  Computing 

4.  Sigal  Gottlieb,  David  Ketcheson  and  Chi- Wang  Shu  Strong  Stability  Preserving  Runge  Kutta 
and  Multistep  Time  Discretizations.  World  Scientific  Press.  January  2011.  ISBN  978-981- 
4289-26-9 

Honors  &  Awards  Received: 

Sigal  Gottlieb  was  awarded  a  chancellor  colloquium  series  lecture  at  UMass  Dartmouth,  which 
included  a  $1,000  award. 


AFRL  Point  of  Contact:  Vassilios  Kovanis,  PhD,  Electro-Optics  Components  Branch 
Sensors  Directorate.  Wright  Patterson  AFB  OH.  Phone  937-904-9943.  Initial  contact  by  email  on 
3/24/09,  talked  twice  on  the  phone  about  the  use  of  SSP  time  discretization.  Arranged  to  meet  in 
Fall  2009. 

Transitions:  Information  about  these  numerical  methods  was  requested  by 
Frank  Giraldo,  (email:  fxgirald@nps.edu,  Assoc  Professor  and  Assoc  Chair  for  Research,  Depart¬ 
ment  of  Applied  Mathematics,  Monterey,  CA)  on  8/12/08. 

Uri  Shumlak,  (Professor,  Aerospace  &;  Energetics  Research  Program  University  of  Washington, 
Seattle,  WA),  who  thought  that  SSP  time-discretization  may  be  applicable  to  his  research  devel¬ 
oping  a  DG  method  for  the  two-fluid  plasma  model.  Sent  him  this  information  on  8/18/08. 
Marsha  Berger,  (Computer  Science  Department,  Courant  Institute  of  Mathematical  Sciences,  New 
York  University)  on  8/14/08. 

Zhilin  Li,  Professor,  (e-mail:  zhilin@math.ncsu.edu,  Department  of  Mathematics,  North  Carolina 
State  University,  Raleigh,  NC)  used  one  of  these  implicit  SSPRK  schemes  for  free  boundary /moving 
interface  problems.  We  discussed  this  by  email  on  11/1/09. 
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